Below is a map of culverts identified in the Washington Department of Fish and Wildlife’s Fish Passage Barrier Inventory as of January 2021. Each barrier has been assigned to a “Cost Percentile” representing the predicted cost of correcting the barrier relative to others in the database. Predictions are based boosted regression trees fit on 1,235 barrier correction projects documented in the Pacific Northwest Salmonid Habitat Projects database. Detailed methods are available in a manuscript currently under peer review and available upon request.

pal <- 
  colorQuantile(
    palette = "RdBu",
    domain = sf_inv_preds$costpred_brt,
    n = 10,
    reverse = TRUE
  )

leaflet(sf_inv_preds, width = "100%") %>%
  addTiles %>% # Add default OpenStreetMap map tiles
  addCircleMarkers(
    radius = 5,
    weight = 1.5, 
    color = "black", 
    opacity = 1, 
    fillOpacity = 1, 
    fillColor = ~ pal(costpred_brt), 
    clusterOptions = markerClusterOptions(spiderfyOnMaxZoom = FALSE, disableClusteringAtZoom = 10),
    popup = 
      ~ paste0(
        "<b>WDFW Site ID:</b> ", 
        "<a href='http://apps.wdfw.wa.gov/fishpassagephotos/Reports/", site_id, "_Report.pdf' target = '_blank'>",
        site_id,
        "</a><br>", 
        "<b>Cost Percentile:</b> ", match(pal(costpred_brt), rev(brewer.pal(10, "RdBu")))*10, "%",
        "<br>",
        "<b>Passability:</b> ", case_when(percent_fish_passable_code == 99 ~ "Unknown", 
                                              percent_fish_passable_code == 10 ~ "0%", 
                                              percent_fish_passable_code == 20 ~ "33%", 
                                              percent_fish_passable_code == 30 ~ "66%", 
                                          ),
        "<br>",
        "<b>Potential Species:</b>", case_when(is.na(potential_species) ~ " ",
                                               TRUE ~ ""), 
        potential_species,
        "<br>",
        "<b>Lineal Gain:</b> ", scales::comma(lineal_gain_measurement, 1), case_when(is.na(lineal_gain_measurement) ~ "",
                                                                                     TRUE ~ "m"),
        "<br>",
        "<b>Ownership:</b> ", case_when(owner_type_code == 1 ~ "City", 
                                        owner_type_code == 2 ~ "County", 
                                        owner_type_code == 3 ~ "Federal", 
                                        owner_type_code == 4 ~ "Private",
                                        owner_type_code == 5 ~ "State", 
                                        owner_type_code == 6 ~ "Tribal", 
                                        owner_type_code == 7 ~ "Other", 
                                        owner_type_code == 8 ~ "Port",
                                        owner_type_code == 9 ~ "Drainage District", 
                                        owner_type_code == 11 ~ "Irrigation District", 
                                        owner_type_code == 12 ~ "Unknown", 
                                        # 1 = "city", 2 = "county", 3 = "federal", 4 = "private", 5 = "state", 6 = "tribal", 7 = "other", 8 = "port", 9 = "drainage district, 11 = "irrigation district", 12 = "unknown"
        ),
        "<br>",
        "<b>Survey Date:</b> ", survey_date,
        ""
      )
  ) %>%
  addLegend(colors = rev(brewer.pal(10, "RdBu")), labels = paste0(c(1:10)*10, "%"), opacity = 1, title = "Cost Percentile") %>%
  # addPopups %>%
  setView(lng = -122, lat = 47.5, zoom = 7)
# Present regional and ownership summaries ----
# __Plot densities ----

sf_inv_preds %>% 
  st_drop_geometry() %>% 
  # filter(Ownership == "County") %>% 
  mutate(
    value = "county"
  ) %>%
  select(
    county_name,
    value,
    costpred_brt
  ) %>%
  bind_rows(
    sf_inv_preds %>% 
      st_drop_geometry() %>% 
      select(
        county_name,
        # value,
        costpred_brt
      ) %>%
      expand_grid(county_name2 = unique(sf_inv_preds$county_name)) %>%
      mutate(
        value = "all",
        county_name = county_name2,
        county_name2 = NULL
      )
  ) %>%
  filter(county_name != "") %>%
  group_by(county_name, value) %>%
  mutate(
    fill_col = case_when(
      value == "all" ~ NA_real_,
      value == "county" ~ mean(exp(costpred_brt))
    ),
    fill_val = mean(exp(costpred_brt))
  ) %>%
  ungroup() %>%
  arrange(-fill_col) %>%
  mutate(
    county_name = fct_inorder(factor(county_name))
  ) %>%
  ggplot() + 
  geom_density(
    aes(
      x = exp(costpred_brt),
      # y = after_stat(scaled),
      fill = fill_col,
      group = value,
      # alpha = value
    ),
    # stat = "scaled",
    alpha = 0.5
    # bins = 100
    # position = position_stack()
  ) + 
  geom_vline(
    aes(
      xintercept = fill_val,
      group = value,
      # color = fill_col
      linetype = value
    ),
    # size = 1,
    # linetype = "dashed",
    # data = ~ .x %>%
    #   group_by(county_name, value) %>%
    #   summarise(mean_cost = mean(exp(costpred_brt))) %>%
    #   filter(county_name != "")
  ) +
  geom_label(
    aes(
      x = 7e3,
      y = 2.5, 
      label = n
    ),
    hjust = 0.05,
    size = 2.75,
    data = ~ .x %>% filter(value == "county") %>% group_by(county_name) %>% summarize(n = paste("N =", scales::comma(n())))
  ) +
  scale_x_log10(
    expression("" %<-% "Lower Cost                 (Log - Scale)                 Higher Cost" %->% ""),    
    breaks = NULL
  ) +
  scale_y_continuous(
    "Density",
    labels = NULL,
  ) +
  # scale_alpha_manual(
  #   values =
  #     c(
  #       "county" = 0.5,
  #       "all" = 0
  #     )
  # ) +
  # scale_fill_manual(
  #   values = 
  #     c(
  #       "county" = "grey",
  #       "all" = "white"
  #     )
  # ) +
  # scale_color_distiller(palette = "RdBu", na.value = "black") +
  scale_linetype_manual(
    values =
      c(
        "county" = "solid",
        "all" = "dashed"
      )
  ) +
  scale_fill_distiller(palette = "RdBu", na.value = "grey90") +
  ggthemes::theme_clean() +
  theme(
    legend.position = "none",
    plot.background = element_rect(color = "white"),
    # panel.backgroun = element_rect(fill = "grey85")
  ) +
  facet_wrap(
    "county_name", 
    # scales = "free_y",
    ncol = 3
  ) +
  labs(
    title = "Predicted cost distribtuions by county",
    subtitle = str_wrap("Solid line represents county mean, dashed line represents overall mean. Background curve represents overall distribution across all counties. Area under all curves scaled to one (i.e., sizes of curves do not represent relative number of barriers).")
  )

# Codes needed (1 = "city", 2 = "county", 3 = "federal", 4 = "private", 5 = "state", 6 = "tribal", 7 = "other", 8 = "port", 9 = "drainage district, 11 = "irrigation district", 12 = "unknown")
sf_inv_preds %>% 
  st_drop_geometry() %>% 
  # filter(Ownership == "County") %>% 
  mutate(
    value = "owner_type_code",
    owner_type_code =
      case_when(
        owner_type_code == 1 ~ "City",
        owner_type_code == 2 ~ "County", 
        owner_type_code == 3 ~ "Federal", 
        owner_type_code == 4 ~ "Private",
        owner_type_code == 5 ~ "State",
        owner_type_code == 6 ~ "Tribal",
        owner_type_code == 7 ~ "Other",
        owner_type_code == 8 ~ "Port",
        owner_type_code == 9 ~ "Drainage district",
        owner_type_code == 11 ~ "Irrigation district",
        owner_type_code == 12 ~ "Unknown"
      )
  ) %>%
  select(
    owner_type_code,
    value,
    costpred_brt
  ) %>%
  bind_rows(
    sf_inv_preds %>% 
      st_drop_geometry() %>% 
      select(
        owner_type_code,
        # value,
        costpred_brt
      ) %>%
      expand_grid(owner_code2 = unique(sf_inv_preds$owner_type_code)) %>%
      mutate(
        value = "all",
        owner_type_code = as.character(owner_code2),
        owner_code2 = NULL
      )
  ) %>%
  filter(!(owner_type_code %in% c(1:12))) %>% drop_na() %>%
  group_by(owner_type_code, value) %>%
  mutate(
    fill_col = case_when(
      value == "all" ~ NA_real_,
      value == "owner_type_code" ~ mean(exp(costpred_brt))
    ),
    fill_val = mean(exp(costpred_brt))
  ) %>%
  ungroup() %>%
  arrange(-fill_col) %>%
  mutate(
    owner_type_code = fct_inorder(factor(owner_type_code))
  ) %>%
  ggplot() + 
  geom_density(
    aes(
      x = exp(costpred_brt),
      # y = after_stat(scaled),
      fill = fill_col,
      group = value,
      # alpha = value
    ),
    # stat = "scaled",
    alpha = 0.5
    # bins = 100
    # position = position_stack()
  ) + 
  geom_vline(
    aes(
      xintercept = fill_val,
      group = value,
      # color = fill_col
      linetype = value
    ),
    # size = 1,
    # linetype = "dashed",
    # data = ~ .x %>%
    #   group_by(county_name, value) %>%
    #   summarise(mean_cost = mean(exp(costpred_brt))) %>%
    #   filter(county_name != "")
  ) +
  geom_label(
    aes(
      x = 7e3,
      y = 2.5, 
      label = n
    ),
    hjust = 0.05,
    size = 2.75,
    data = ~ .x %>% filter(value == "owner_type_code") %>% group_by(owner_type_code) %>% summarize(n = paste("N =", scales::comma(n())))
  ) +
  scale_x_log10(
    expression("" %<-% "Lower Cost                 (Log - Scale)                 Higher Cost" %->% ""),    
    breaks = NULL
  ) +
  scale_y_continuous(
    "Density",
    labels = NULL,
  ) +
  # scale_alpha_manual(
  #   values =
  #     c(
  #       "county" = 0.5,
  #       "all" = 0
  #     )
  # ) +
  # scale_fill_manual(
  #   values = 
  #     c(
  #       "county" = "grey",
#       "all" = "white"
#     )
# ) +
# scale_color_distiller(palette = "RdBu", na.value = "black") +
scale_linetype_manual(
  values =
    c(
      "owner_type_code" = "solid",
      "all" = "dashed"
    )
) +
  scale_fill_distiller(palette = "RdBu", na.value = "grey90") +
  ggthemes::theme_clean() +
  theme(
    legend.position = "none",
    plot.background = element_rect(color = "white"),
    # panel.backgroun = element_rect(fill = "grey85")
  ) +
  facet_wrap(
    "owner_type_code", 
    # scales = "free_y",
    ncol = 3
  ) +
  labs(
    title = "Predicted cost distribtuions by ownership",
    subtitle = str_wrap("Solid line represents ownership class mean. Area under all curves scaled to one (i.e., sizes of curves do not represent relative number of barriers).")
  )